function [Teci2crtbp,Tcrtbp2eci,S_eci] = CRTBP2ECI(JD_interest,S_crtbp,JD_EM, Rme, Vme)

% GM_sun=1.327124400*10^(11);     %Sun
GM_earth=3.986004418*10^(5);    %Earth
GM_moon = 4.9048695*10^(3);     %Moon
%[km^3/s^2] Gravitational Parameters

r12 = 384400;
%[km]Distance between the Earth and the Moon.

Wcm = sqrt((GM_earth+GM_moon)/r12^3)*[0;0;1];
%[] angular velocity in Cr3bp frame

TU = 3.75676967e5;
LU = r12;
VU = LU/TU;
%[] Canonical unit conversions

S_crtbp = [S_crtbp(1:3)*LU;S_crtbp(4:6)*VU];
R_crtbp = S_crtbp(1:3);
V_crtbp = S_crtbp(4:6);
%[] Non-inertial position and velocity of the spacecraft at JD_interest.
%% 1: Using julian date, determine moon's state at interested julian date.

JD_difference_vector = JD_EM-JD_interest;
%[] Determine the difference time vector

ind = find(JD_difference_vector > 0 , 1 , 'first') - 1;
%[] Find the index of most close previous state

Rmoone_o = Rme(:,ind);
Vmoone_o = Vme(:,ind);
Smoone_o = [Rmoone_o;Vmoone_o];
%[km,km/s] Position and velocity of the moon WRT Earth.

UTC_JD_interest = CalendarDate(JD_interest);
UTC_JD_closest = CalendarDate(JD_EM(ind));
%[] Get UTC time of JD time of interest, which is the JD_interest and
%JD(ind), which is the closest previous information.

% time_difference_seconds = etime(UTC_JD_closest,UTC_JD_interest);
time_difference_seconds = etime(UTC_JD_interest,UTC_JD_closest);

%[] Determine time difference in seconds.
if time_difference_seconds == 0
    Rmoone_interest = Rmoone_o;
    Vmoone_interest = Vmoone_o;
    Smoone_interest = [Rmoone_interest;Vmoone_interest];
else
    [~,S] = ode45(@(t,S)R2bpCartesianModel(t,S,GM_earth),[0,time_difference_seconds],Smoone_o);
    S = S';
    %[] Integrate for difference amount to determine exact state of the moon
    %wrt Earth in ECI frame.
    
    Rmoone_interest = S(1:3,end);
    Vmoone_interest = S(4:6,end);
    Smoone_interest = [Rmoone_interest;Vmoone_interest];
    %[] Position and velocity of the moon at JD_interest.
end
%% 2: Determine states of CM

Scme = (GM_moon*Smoone_interest)/(GM_moon+GM_earth);
Rcme = Scme(1:3);
% Vcme = Scme(4:6);
%[] Determine the state of the center of mass at JD of interest
%There are no additional terms because the Earth's position wrt Earth is
%zero vector and the gravitational parameter of the satellite is assumed to
%be zero.

%% 3: Make all states WRT CM

Searthcm = -Scme;
% Rearthcm = Searthcm(1:3);
% Vearthcm = Searthcm(4:6);
%[] State of the Earth wrt CM

Smooncm = Smoone_interest-Scme;
% Rmooncm = Smooncm(1:3);
% Vmooncm = Smooncm(4:6);
%[] State of the moon wrt CM

% Ssatcm = S_eci-Scme;
% Rsatcm = Ssatcm(1:3);
% Vsatcm = Ssatcm(4:6);
%[] State of the state of interest [S_eci] wrt Center of mass.

%% 4: Determine Teci2crtbp and Tcrtbp2eci

Hm = cross(Rmoone_interest,Vmoone_interest);
%[1/s]Current specific angular momentum of the Moon WRT the CM in ECI coordinates.

Rhat = Rmoone_interest / norm(Rmoone_interest);
%[]Current CRTBP radial direction in ECI coordinates.

Nhat = Hm / norm(Hm);
%[]Current CRTBP normal direction in ECI coordinates.

That = cross(Nhat,Rhat);
%[]Current CRTBP tangential direction in ECI coordinates.

Tcrtbp2eci = [Rhat, That, Nhat];
%[]Matrix that transforms vectors from CRTBP to ECI coordinates.

Teci2crtbp = transpose(Tcrtbp2eci);
%[] Matrix that transforms vectors from ECI to CRTBP coordinates.


%% 5: Determine S_crtbp

Rsatcm_eci = Tcrtbp2eci * R_crtbp;
Vsatcm_eci = Tcrtbp2eci*(V_crtbp+cross(Wcm,Teci2crtbp*Rcme)+cross(Wcm,R_crtbp));
Ssatcm_eci = [Rsatcm_eci;Vsatcm_eci];
S_eci = Ssatcm_eci + Scme;

end